Hmmrank coordinates - #72
Conversation
Both pinned to nf-core/modules 106170c99, which adds alignment coordinates to the ranked table: hmmsearch's per-domain output is now passed on to the ranking step, and when it is present each row gains the sequence and profile lengths plus bounds, coverage and island count for the hmm, ali and env coordinate sets. Only the ranked table's checksum moves in the pipeline snapshot. The per-profile fastas SEQTK_SUBSEQ extracts are byte-identical, so which sequences are selected is unaffected -- the new columns are additive, and significance is still decided by the full-sequence scores. Generated by Claude
Describe the columns in docs/output.md, and point at them from the --save_domtblout section of docs/usage.md, since the ranked table is an easier place to read coordinates off than the per-domain table itself. Also correct two claims in AGENTS.md that had HMMER_HMMRANK and both fasta_* subworkflows down as pipeline-local. All three are vendored nf-core/modules components tracked in modules.json; they originated here and are still maintained by this pipeline's author, which is what made the wording plausible. Since nf-core/ampliseq vendors hmmer/hmmrank and fasta_hmmsearch_rank_fastas too, note that changes to them have to stay backwards compatible. Generated by Claude
This PR is against the
|
Generated by Claude
|
❌ nf-test failed with latest Nextflow versionNote Tests with Nextflow's latest version failed but it will not cause a CI workflow failure.
See the full run for details. |
piplus2
left a comment
piplus2
left a comment
There was a problem hiding this comment.
I've asked Claude to look at the modifs:
Findings
1. Duplicate rank rows when an HMM file holds multiple models
Severity: medium
Location: modules/nf-core/hmmer/hmmrank/main.nf:68
domain_coords <- domains %>% distinct(accno, profile, tlen, qlen)This assumes one row per (accno, profile), but profile is derived from the file name, and a single *.domtbl.gz can hold hits from several query profiles when the supplied HMM file contains more than one model — which this pipeline permits, since extract_hmm in the phylosearch samplesheet is optional and without it the whole HMM file goes to hmmsearch as-is.
Different models have different qlen, so distinct() emits one row per model for the same (accno, profile) key. Confirmed in R: with two qlen values for one key, distinct() returns 2 rows.
Impact. The left_join(domain_coords, by = c('accno','profile')) at lines 74/98 then multiplies every matching rank row (N tblout rows × M models), so *.hmmrank.tsv.gz gains duplicate rows. The subworkflow's filter { it.rank == '1' } → collectFile writes the same accession several times into the SEQTK_SUBSEQ id file, duplicating extracted query sequences that go on into placement. Before this PR the same input produced N rows, so this is a regression.
islands() also silently unions domains belonging to different models into one coordinate range.
Fix. Aggregate domain_coords to exactly one row per key, or join on profile_desc (the query name) as well — either closes both problems.
2. NA coordinate columns for tblout-only hits
Severity: low
Location: modules/nf-core/hmmer/hmmrank/main.nf:74
The left_join leaves all 14 new columns as NA for any (accno, profile) present in the tblout but absent from the domtblout. This is reachable with stock thresholds: hmmsearch reports a sequence in --tblout on the sequence reporting threshold, while --domtblout rows require the domain reporting threshold (--domE, default 10). A hit whose sequence score clears but whose every individual domain is above --domE yields zero domain rows.
docs/output.md and meta.yml both state each row carries the columns, with no mention of NA, and the new nf-test assertion ...isInteger() on hmm_len would fail on such a row rather than reporting it.
3. Eager output extraction masks process failures
Severity: low
Location: modules/nf-core/hmmer/hmmrank/tests/main.nf.test:78 (and :145)
def rows = path(process.out.hmmrank[0][1]).linesGzipThis is evaluated eagerly, before assert process.success inside assertAll. If the process fails, process.out.hmmrank is empty and this throws IndexOutOfBoundsException from the then block, masking the actual process failure and its log.
Fix. Move the extraction inside the assertion closures, or assert success first.
4. rows[1..-1] silently misbehaves on header-only output
Severity: low
Location: modules/nf-core/hmmer/hmmrank/tests/main.nf.test:147
rows[1..-1] on a header-only output does not throw: in Groovy, 1..-1 on a 1-element list resolves to the reverse range 1..0 and returns the header row itself, so the "known table" comparison runs against a bogus row list instead of failing cleanly. The sibling test guards this with assert rows.size() > 1 (line 90); this one does not.
Checked and correct
- The
separate()field mapping against HMMER's 23-column domtblout layout - The
+ 1Labutment rule andlag(cummax(t), default = 0L)island detection - Stub header column order matching the real
left_joinoutput order - The
^prefix\./\.tbl\.gz$anchoring — a strict improvement, since the unanchored.tbl.gzcould match inside a profile name; the no-domtblout snapshot md5 is unchanged - The AGENTS.md local-vs-vendored claim:
modules/local/hmmerandsubworkflows/local/utils_nfcore_phyloplace_pipelineare indeed the only local components
Verification performed
The island/union arithmetic was checked in R against synthetic inputs, and the new channel plumbing (collect/ifEmpty/combine) in an actual Nextflow run — including the empty-domtblout case producing [meta, tblouts, []] and the Groovy-falsy branch in the script.
A hit is reported on the per-sequence threshold while a domain has to clear the per-domain one, so a ranked row can have no domtblout records behind it. Keeping NA there is intentional -- it is the clearest signal that hmmsearch reported no domains for the hit -- so say so in the docs. Raised in review of nf-core#72. The remaining review findings are upstream in the vendored module and are tracked in nf-core/modules#12787. Generated by Claude Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PDswY3eJxbJ9xnWSCQFVRG
|
Thanks @piplus2 — useful review. I went through the findings with some help from Claude Code, and here's where each one landed. The important structural point first: all four findings are in Finding 1 — real, and worth fixing, but a corner case. Confirmed the preconditions: Finding 2 — correct reading, but the Finding 3 — real, minor. Agreed, lines 78 and 145 read the output before Finding 4 — doesn't hold. I checked the Groovy behaviour rather than reasoning about it: Your "checked and correct" list matches what I'd verified independently, including the |
|
Thanks @piplus2. Since your main point -- an important finding -- is meant to be fixed in modules, would you care to look at nf-core/modules#12788? |
|
Ok checking in a few! |
qlen is the profile's own length, so hmm_len / qlen is the proportion of the profile a hit covers and the hmm bounds say how much is missing from each end. Everything needed is already on the row; the columns just did not say what they were for. Generated by Claude Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PDswY3eJxbJ9xnWSCQFVRG
tlen is the target sequence's length, so ali_len / tlen is the proportion of the sequence the alignment takes up, the counterpart to hmm_len / qlen for the profile. Spell out which length goes with which coordinate set, since hmm positions live in the profile while ali and env positions live in the sequence, and dividing across the two spaces means nothing. Generated by Claude Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PDswY3eJxbJ9xnWSCQFVRG
Pulls in nf-core/modules#12788. An HMM file may hold several models, in which case one search puts them all in one table under one file-derived profile name, and the coordinate columns keyed on that alone: rows were multiplied and coordinates unioned across different models. Reachable here because extract_hmm is optional, so a multi-model file goes to hmmsearch as-is. Also breaks rank ties on the model name, which matters because rank 1 is what the subworkflow selects sequences on. Pipeline snapshot unchanged: all test HMMs that skip extract_hmm hold a single model, and the fix is a no-op in that case. Generated by Claude Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PDswY3eJxbJ9xnWSCQFVRG
* hmmer/hmmrank: key domain coordinates on the query model The profile column comes from the file name, but an HMM file may hold several models, and the whole file is then searched in one go, so one table can carry rows for several query names under one profile. The coordinate columns added in nf-core#12781 keyed on (accno, profile) only, which that breaks in two ways: distinct() returns one row per model, so the left_join against a tblout side with one row per model becomes a cartesian product, and islands() unions domains across different models, which is meaningless for hmm_from/hmm_to since those are positions in the profile. Concatenating the two existing NrdJ fixtures into one file reproduces it: 17 rows instead of 9, and every coordinate corrupted -- NrdJ's hit on AACBMDAH_04658 reports NrdJm1's 1-397 rather than its own 11-95, and AAGAFPGM_00302 gets 6-711 with hmm_len 676, a cross-model union matching neither model. The query name is already on both sides -- column 3 of the tblout, parsed as profile_desc, and column 4 of the domtblout, until now discarded as a dummy -- so carrying it through islands() and the join fixes both halves. Output columns are unchanged: the join maps profile_desc onto the domtblout query, so no new column is emitted and single-model input produces byte-identical output. All existing module and fasta_hmmsearch_rank_fastas snapshots pass untouched. Adds a regression test that builds the multi-model case by gunzipping, concatenating and re-gzipping the two fixtures inline, so no near duplicate test data is needed. It fails on the previous code. Also from the review of nf-core/phyloplace#72, both in tests/main.nf.test: output was read in the then block before assert process.success ran, so a process failure surfaced as an IndexOutOfBoundsException with its log hidden; and the hmm_len assertion required an integer on every row, ruling out the NA that a hit with no domain records legitimately gets. NA is intended -- it is the clearest signal that hmmsearch reported no domains for a hit -- so the assertion now allows it and meta.yml says so. Closes nf-core#12787 Generated by Claude Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PDswY3eJxbJ9xnWSCQFVRG * hmmer/hmmrank: break rank ties on the model name too Review follow-ups on nf-core#12788. The rank tie-break had the same root cause as the join: arrange() fell back on profile, which is the file name, so for a file holding several models it was constant across them and ties silently resolved to input order. rank 1 is what downstream consumers select on, so that made the selection arbitrary. Adding profile_desc restores the documented alphabetical tie-break; it only takes effect where score, evalue and profile are all equal, which cannot happen across separate files, so single-model input is unaffected and no snapshot moves. The hmm_len assertion allowed NA on every row, which is exactly what a broken join produces, so it had stopped guarding the regression closest to it. It now also requires at least one integer. meta.yml said one row per (accno, profile) pair, stale after this PR, and profile_desc becomes load-bearing without its name saying that it holds the query model name. Both spelled out now. Also close the input streams in the fixture-combining closure. Generated by Claude Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PDswY3eJxbJ9xnWSCQFVRG --------- Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
2 participants
PR checklist
nf-core pipelines lint).nextflow run . -profile test,docker --outdir <OUTDIR>).nextflow run . -profile debug,test,docker --outdir <OUTDIR>).docs/usage.mdis updated.docs/output.mdis updated.CHANGELOG.mdis updated.README.mdis updated (including new tool citations and authors/contributors).Description
Closes #70.
Adds alignment coordinates and lengths to the ranked
hmmsearchsummary, so a hit can be judged on how much of the profile it actually covers rather than on score alone.Third and last PR of the chain: nf-core/test-datasets#2232 added the fixtures, nf-core/modules#12781 did the work in
hmmer/hmmrankandfasta_hmmsearch_rank_fastas,and this pulls both into the pipeline and documents them.
What changes
When
--save_domtbloutis set,*.hmmrank.tsv.gzgains 14 columns:tlenandqlen, plusfrom,to,lenandn_islandsfor each of thehmm,aliandenvcoordinate sets.lenis the size of the union of that hit's domains, so it drops belowto - fromwhen the domains are scattered, andn_islandscounts how many separate stretches they form.Coverage spans every domain hmmsearch reported, however weak on its own. That is deliberate: a sequence can clear the full-sequence threshold while each of its domains is poor in isolation, and those fragmented hits are exactly the ones these columns exist to expose. Which sequences are reported at all is still decided by the full-sequence scores, unchanged from before.
Without
--save_domtbloutthe table is exactly as it was.Also here
Two claims in
AGENTS.mdhadHMMER_HMMRANKand bothfasta_*subworkflows down as pipeline-local. All three are vendored nf-core/modules components tracked inmodules.json— they originated here and are still maintained by this pipeline's author, which is what made the wording plausible, but they must be changed upstream rather than edited in place. Corrected, with a note that nf-core/ampliseq vendors two of them so changes have to stay backwards compatible.Testing
nf-test6/6,nf-core pipelines lint0 failures,prek run -aclean,nextflow lint .clean at 25.10.4 and 26.04.6 (remaining findings are pre-existing ondev).Only the ranked table's checksum moves in the pipeline snapshot — the per-profile fastas
SEQTK_SUBSEQextracts are byte-identical, confirming the new columns don't affect selection.